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Abstract 


This report investigates two of the most common modes of localized failures; namely, 
periodic fiber-bridged matrix cracks, and transverse matrix cracks. A modification of 
Daniels’ bundle theory is combined with Weibull’s weakest link theory to model the 
statistical distribution of the periodic matrix cracking strength for an individual layer. 
Results of the model predictions are compared with experimental data from the open 
literature. Extensions to the model are made to account for possible imperfections wit hin 
the layer (i.e., non-uniform fiber lengths, irregular crack spacing, and degraded in-situ 
fiber properties), and the results of these studies are presented. A generalized shear-lag 
analysis is derived which is capable of modeling the development of transverse matrix 
cracks in material systems having a general multilayer configuration and under states of 
full in-plane load. A method for computing the effective elastic properties for the damaged 
layer at the global level is detailed based upon the solution for the effects of the damage 
at the local level. This methodology is general in nature and is therefore also applicable 
^ systems. The characteristic stress-strain response for more general cases is 

shown to be qualitatively correct (experimental data is not available for a quantitative 
evaluation), and the damage evolution is recorded in terms of the matrix crack density as 
a function of the applied strain. Probabilistic effects are introduced to account for the 
statistical nature of the material strengths, thus allowing cumulative distribution curves for 
the probability of failure to be generated for each of the example laminates. Additionally, 
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Abstract 


Oh and Finney’s classic work on fracture location in brittle mat^ials is extended and 
combined with the shear-lag analysis. The result is an analytical form for predicting the 
probability density function for the location of the next transverse crack occurrence within 
a crack bounded region. The results of this study verified qualitatively the validity of 
assuming a uniform crack spacing (as was done in the shear-lag model). 
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Chapter 1 


Introduction 


Laminated composite materials can be engineered to exhibit gradual (or progressive) 
failure as opposed to the less desirable catastrophic failure. Localized failures begin to 
occur early in the loading history of a laminate, but due to the laminate’s ability to 
redistribute its load internally, the laminate as a whole remains unfailed and continues to 
perform its load-carrying function. As the loading of the laminate continues, the number 
of localized failures accumulate and begin to degrade the overall laminate properties, 
eventually reaching the point at which total laminate failure occurs. 

In trying to model laminate failure, it is important to recognize the progressive nature 
of the failure process in order to obtain a correct strength prediction. Ignoring this 
behavior, as is done in first ply failure (FPF) theories, leads to failure predictions which 
are too conservative. When implementing a progressive failure scheme, the unit of failure 
which is to be considered must be defined. At the macroscopic level, this unit can be 
assumed to either represent the entire load bearing capacity of a ply ("ply failure"); or this 
ply can be further sub-divided into smaller units, each corresponding to a particular load 
bearing mode ("modal failure"). Techniques allowing for such behavior have been pursued 
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as an improvement over FPF (Petit and Waddoups, 1969; Thomas and Wetherhold, 
1991(b)). HowevCT, these methods also are consCTvative because they do not allow for any 
load carrying capacity of the damaged layer away from the localized failure zone. The 
load redistribution scheme utilized in this iqmrt employs the modal approach with the 
belief that it better represents the progressive frilure nature of the composite. A non- 
interactive failure function is used in which the individual stresses are assumed to act 
indq>endently of one another in the failure process. Four potential modes of failure are 
recognized within each layer (see Figure 1.1). In the longitudinal direction a two-stage 
failure process is considered. Both of these modes of failure are associated with the in- 
plane normal stress in the fiber direction, a ^ . The first stage is representative of the 
periodic matrix cracking (cracks running perpendicular to the fiber direction and bridged 
by the fib^) predicted by ACK theory. Failure of this mode degrades the stiffness 
properties of the layer, but the survival of the fibers which bridge these cracks allows 
continued loading. The second stage of failure deals with the fracture of these bridging 
fibers. In the transverse direction, a single mode of failure models the occurrence of 
transverse matrix cracking (through-the-thickness cracks running parallel to the fiber 
direction). This failure mode is a function of the normal stress, . The fourth and final 
mode is an in-plane shear failure associated with the stress Xjj (o^ in Voigt-contracted 
notation). The failure functions for each of these modes are given below in equations (1.1) 
through (1.4). 
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Periodic, bridged matrix cracks 


Fiber failure 



Transverse matrix craddng 



Figure 1.1: Potential failure modes for an individual layer. 
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where 



( 1 . 2 ) 



(1.3) 



(1.4) 


<*) 


X 



0 , xiO 
X , x>0 


(1.5) 


and Xj and xf are the tensile and compressive strengths, respectively. Values of / less 
than one (/ < 1) designate the safe loading regime, and values of / equal to one(/ - 1) 
denote failure. (Keep in mind that the shear strengths must be equal, X^ - X^)- 

Unlike conventional engineering materials metals), where local failures can be 
rdieved by plastic deformation thus allowing their strragths to be considered determimstic- 
ally, composite materials’ brittle nature results in thdr being characterized by a high 
variability of material strengths. Therefore, in order for effidait appUcation of these 
materials, design methods using probabilistic fmlure analyses must be advanced. As a 
result, rather than having a clearly defined fail/m fail situation as described in the 
previous paragraph, the situation becomes a matter of evaluating the probability of 
failure/survival ( i.e., Pr{f <: 1 )) for a given system where the material’s strengths are 
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considered to be random variables. Within this report, investigations will be made into 
modeling both the progressive and probabilistic aspects of failure. 

With regard to failure in the form of periodic matrix cracking, the use of a modified 
bundle theory to model the post-matrix cracking state predicted by ACK theory has been 
proposed (Walls, 1986; Evans and Marshall, 1989). Chapter 2 will present a scheme for 
incorporating this method into a laminate analysis algorithm which utilizes Monte Carlo 
simulations to evaluate reliability. Example calculations using this method will be 
demonstrated and compared with experimental results from the open literature. Some of 
the problems encountered with the method will be discussed and potential causes will be 
investigated. 

The onset of transverse matrix cracking has been foimd to be a key occurrence in the 
laminatp- fail ure process, and as a result much research has been devoted to its modeling. 
Continuum damage models (Talreja, 1985; Nuismer and Tan, 1988) as well as many shear- 
lag approaches (Garret and Bailey, 1977; Reifsnider, 1978; Laws and Dvorak, 1988; Lee 
and Daniel, 1990) have been proposed for analyzing cross-ply laminates. Chapter 3 
extends the shear-lag method of Lee and Daniel from a two-layer cross-ply system to a 
general symmetric multilayer system. In Chapter 4, the elasticity problem for the region 
of the laminate between two parallel matrix cracks having a general off-axis orientation is 
posed from equilibrium considerations. This is done in terms of the average (through-the- 
thickness) stresses and solved using the generalized shear-lag relation and the appropriate 
boundary conditions. A method for modeling the effective elastic behavior of the damaged 
material based on this solution is detailed. This model is then included in a computer 
program designed for probabilistic laminate analysis and the results are compared to those 
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determined using the ply drop-off technique. 

In Chapter 5, a cross-ply laminate experiencing transverse matrix cracking under 
uniaxial longitudinal loading is studied. Combining the solution for the stress state in the 
damaged layer (determined in Chapter 4) with Oh and Finnie’s results for fracture location 
in brittle solids (Oh and Finnie, 1970) an expression is determined for the probability 
doisity function for the location of the occurrence of the next transverse matrix crack. 
Example problems are presented, and the validity of assuming a uniform crack spacing in 
models such as that of Chapter 4 is examined qualitatively. 
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Periodic, Bridged 
Matrix Cracks 


The use of a modified bundle theory to model the post-matrix cracking state predicted 
by ACK theory has been proposed (Walls, 1986; Evans and Marshall , 1989), In this 
chapter the suitability of this method is examined by incorporating it into a load sharing 
algorithm for the progressive failure study of composite laminates. Example problems will 
be presented where both the modified bundle approach as well as a previous version of the 
load sharing algorithm (Thomas and Wetherhold, 1991(b)) which did not permit separate 
considraation of periodic matrix cracking have been used. The results of these anal yses 
will be compared to experimental results from the open literature. Some of the problems 
encountered with the method will be discussed and potential causes will be investigated. 


LAMINATE ANALYSIS 


Consider an n -layer composite laminate with 4it potential failure modes. As a ramp 
loading is slowly applied some of these modes wiQ begin to fail in some layers. As a 
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lesult of these fsilures, internal loads must be redistributed and stifftiesses reduced. In 
order to model these occurrences, load sharing rules must be adopted that are both 
physically plausible and computationally manageable. Failure is assumed to be governed 
by equations (1.1) through (1.4) of Chapter 1. In the following discussion, a local load 
sharing scheme will be applied by which the load previously carried by the now-failed 
mode will be distributed evenly between the two immediately adjacent layers within the 
laminate. If either of these layers has already failed or contain modes which have failed, 
that portion of the load is globally redistributed in accordance with the laminate constitutive 
law (see Jones (1975) for a review composite mechanics). Upon the occurrence of penodic 
matrix cracking in a particular layer, the and terms of the layer’s stiffnesses 
matrix must be reduced. It is proposed thnt this reduction be made by the ratio of 
corresponding strengths for that layer, i.e.. 


reduced' 



J = 1.2 


( 2 . 1 ) 


For all other modes of failure, the ai^rqpriate stifhiesses are conservatively reduced to 
zero. That is to say, in the case of fiba: failure and Qjj are set to zcto, for 
transverse matrix cracking and are set equal to zero, and finally for the case of 
in-plane shear failure <?66 is set to zero. 


STRENGTH EVALUATION 


At this point the ideal approach would be to state the problem in terms of a failure tree. 
However, the number of different branches quickly becomes too large, and any change in 
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the laminate configuration would require the entire tree to be rebuilt. The use of computer 
simulations is a more viable alternative method of solution. A sample population of 
laminates may be created within the computer’s memory using realizations. That is, 
hypothetical sample values of the modal strengths whose frequency of occurrence is 
weighted in accordance with their distribution. The reliability of the laminate can then be 
determined for a specified load via a Monte Carlo analysis. 

The strengths ^ and are random variables assumed to be characterized by a 
two-parameter Weibull distribution, the parameters of which must be determined 
experimentally. By inverting these distribution functions and sampling a uniform random 
variable on [0,1] the needed realizations for these strengths can be obtained. A rimiiar 
method could also be employed for the strength X^^. However, it has been proposed that 
the charactraization of this strength can be arrived at through more physically based 
arguments via the application of a modified bundle theory. 

After matrix cracking has occurred, the longitudinal strength of the ply is governed by 
the intact fibers which bridge the matrix crack sites. Daniels (1945) has characterized the 
strength of a bundle of fibers under a tensUe axial load as having a normal distribution with 

( 2 . 2 ) 

o = Sf)/nb(sy)[l-b(s^)] (2.3) 

where and a are the average ultimate load and standard deviation of load, respectively, 
for the bxmdle; with n being the number of fibers, b(s) the cumulative distribution 
function for the failure of a single fiber under load s , and is the fiber load which 
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maximizes the expression s[l - b(s)]. For the problem at hand, a modification to this 
theory has been proposed (Walls, 1986; Evans and Marshall, 1989) in which the bridging 
fibers are to be analogous to Daniels’ fibor bundle, with the added compl^ity of the 
fibers’ being embedded in a periodically cracked matrix, thus producing an axial variation 

in the fiber stress. 

Assuming the stress to be constant over the cross sectional area of the fiber, A , the 
probability of failure for the fiber may be described by a Wdbull distribution; 

where a is caUed the WeibuU shape parameter and is dimensionless, and P, is a scale 
parameter based on length. The value of P, is not generally measured directiy from 
experiment. Conventionally, tensile failure data is obtained for fiber specimens all having 
the same experimental test volume, V^, and with each q>ecimen under a state of uniform 
uniaxial tensile stress. This data is fit to a WdbuU distribution of the form 

b(a) = 1 - j] J « * P > 0 

The Wdbull scale parameter in equation (2.5) is based on the experimental volume, and 
is inferred directly from the experimental data. From this value of P , the parameter P, 
can be obtained from the expression below. 



The parameter P has dimensions of [stress] and those of P| are [stress] • [length] . 
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1 

▼s 


Figure 2.1: Axial normal stress distribution within fiber. 

For the case of a fiber which bridges a matrix crack (see Figure 2.1), with the shear 
stress (t) at the debonded fiber/matrix interface assumed to be constant, the axial stress 
distribution within the fiber will be linear. The maximum values will occur at the matrix 

crack faces and minimum values midway between the cracks. Evaluating the stress 
through a force balance. 


a{s,x) = 

A r 


0 i JC i ^ 
2 


(2.7) 


Here ^ is defined to be the fiber load at the matrix crack fece, x is the longitudinal 
position measured from the crack face, r is the fiber radius and d is the crack spacing. 
The value of d must be between x' and 2x' (Aveston et al., 1971), where 
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X 



( 2 . 8 ) 


with IVf ,r,) being the (fiber.mattK) volume ftaction and o„ the cntical matrix cracking 

areas. Itwinbeconservaliyelyassumedthmthevalueofd maintains a uniform value of 

X-. Using the above expression (equation (2.7)) for the fiber stress in equation (2.4) leads 
to the following form of the cumulative distribution function for a single fiber. 


P,Lr 

f— 


^’1 

(2.9) 

Td(a+1) 


rP,j 

‘Pj II 



Performing the necessary maximization results in an impUdt expression for s, . equation 
(2.10), which may be solved numerically. 





( 2 . 10 ) 


TOs may in turn be used in equations (2.2) and (2.3) for the mean strength and standard 
deviation to characterize X^, the post-matrixsnacking ply strength. 


examples 


The system for the example problems was a silicon carbide fiba-reinforced / 

reaction-bonded silicon nitride (SiCVRBSN) tested by Bhatt and Phillips (1990). The test 
specimens measured 12.7 mm wide, had a gage length of 25 mm. and laya thicknesses 
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TABLE 2.1 

Weibull Strength Parameters for SiC/RBSN 
(Bhatt and Phillips, 199^ 


Modal strength, X, 

a 

PQdPd) 

^1. 

6.5 

244. 


5.2 

741. 


10.9 

28. 

x; 

7.5 

56. 

« assumed values 



TABLE 2.2 

Elastic Material Properties for SiC/RBSN 
(Bhatt and Phillips, 1990) 


£j 193. GPa 

69. GPa 

Gu 31. GPa 


of 0.25 mm. The fiber volume fraction was 0.3 and the fiber radius was 71 x 10^ m. The 
Weibull strength parameters are printed in Table 2.1. They were calculated mathematically 
from the published values of the respective means and standard deviations. Table 2.2 
contains information regarding the elastic material properties. Data was also given for the 
individual fibers. Twenty (20) separate tensile tests were performed on individual 
nitrogen-treated SiC fibers at a gage length of 25 mm. The average tensile strength was 
calculated to be 2860 MPa with a standard deviation of 440 MPa and a Weibull shape 
parameter of 8.2. From these results a value of 3043 MPa may be inferred for the 
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WeibuU scale parameter on a unit volume basis, P . 

The results of including the modified bundle model in a laminate analysis are shown in 
Figure 2.2 for a [Oj/POj]^ laminate and in Figure 2.3 for a [+452/-452]^ laminate. In 
addition to these results, each figure also contains results from a similar analysis with the 
excq)tion that periodic matrix cracking was n^lected 0.e., only three potential modes of 
failure per layer, corresponding strengths were )> as wdl as a curve fit to 

e7q>erimentally measured values of the laminate’s mean strength and standard deviation 
(Bhatt and Phillips, 1990). For the [ 02 /^ 2 ]^ laminate, the curve for the modified bundle 
technique shows a very tight distribution and does not agree well with the other two 


14 



Chapter 2: Periodic, Bridged Matrix Cracks 



50000 100000 190000 300000 350000 300000 

Applied Load. N» (N/m) 

Figure 2.3: Comparison of failure mod^s for [+452/-45J, laminate. 


curves. A close examination of the results of the simulation that produced this curve found 
that the mode of failure which ultimately governed the overall failure of the laminate was 
the longitudinal fiber failure of the 0* layers. Furthermore, the modified bundle prediction 
for the coefficient of variation (COV) of was j^roximately 0.05. This value is much 
lower than the COV value of 0.22 observed experimaitally for The combination 

of these two factors led to the steep slope of Figure 2.2. For the [+45^/ ~45j laminate 
the mode of failure which governed the laminate failure was the in-plane shear failure of 
the layers. Thus the low predicted COV of X^ did not effect these results. No 
experimental values for the WeibuU parameters describing the distribution of the shear 
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strraigth Xg were available. The value of the Wdbull scale parameter for was chosen 
to be twice that of the scale paiamet^ for J!^, and a parametric study was done using 
'‘reasonable " values for the WeibuU shape parametm'. The shj^ parameter value that gave 
the best fit, as compared against the ^cperimental data, of the reliability curve for the 
J +45j/ “dSj] laminate was used in Figure 2.3 and is the value listed in Table 2.1. 


NON-IDEAL BUNDLE 


In an effort to e:q)lain the low predicted value for the coefficient of variation of X^ 
obtained using the modified bundle theory, the effects of certain non-ideal conditions which 
may have «dsted, but were not considered in the model, are investigated. While the 
presoit theory has accounted for a variability in the strength from fiber to fiber, it has not 
allowed for any variability in the geometry of the problem. Kerens (1988) has stated that 
all composite fibers are not identical and that these non-uniformities can produce 
unintended but desirable characteristics such as increased toughness. Specifically, Kerans 
studied the effect of non-uniform fiber lengths on the bundle’s stress-strain behavior. To 
find the effect non-uniform fiber lengths would have on the COV of the layer’s ultimate 
stroigth, a computer program was used to simulate a strain controlled test on a bundle of 
such fibers. The number of fibers in the bundle was the same as in the experimental test 
specimens, n =60 . It was assumed that the bundle represented the complete load bearing 
capacity of the damaged layer, and thus the layer’s strragth was calculated by dividing the 
maTimiitn bimdle load by the area of the layer. A series of Monte Carlo analyses was 
p^ormed. A sample population of fiber lengths was simulated for each analysis. Each 
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sample had a normal distribution which was truncated at arbitrarily imposed minimum and 
maximiim values lying five (5) standard deviations away from the mean. For each case 
the minimum value was taken to be the nominal bundle length, and various percentages of 
this value were used for the maximums. The results are shown in Figure (2.4) and in 
Table 2.3. It can be seen that while non-uniform fiber lengths do increase the coefficient 
of variation, the magnitude achieved is still much smaller than that observed experimentally 
and at the same time an adverse effect on the mean strength prediction is experienced. 

Until now it has been assumed that the damaged layer has a regular crack spacing of 
x' . More realistically the crack spacings will be uniformly distributed between the 
tninimiim and maximum values of x' and 2x' allowed by ACK theory. To study what 
effect this would have, simulations were conducted keq)ing the same gage length (L) as 
in the other analyses but allowing the crack spacings (d,) to take on random values 
uniformly distributed between x' and 2x'. That is to say, 

I = ; d,€{x\2x') (2.11) 

i-l 

where n is the number of segments necessary to span the entire gage length. Applying 
the modified bundle theory to a simulation of such a crack field lowered the COV of 
from 0.0453 to 0.0449. Based on this negligible change, it appears that the assumption 
of a uniform crack spacing of x' is reasonable and errs on the consCTvative side. 

Lastly, the effect of a degraded in-situ Weibull shape parameter for the fiber was 
investigated. A parametric study was done over a range of values of from a minimum of 
2.0 to a maximum value of 8.2 which was measured experimentally. The COV ofX^^ 
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TABLE 2.3 

Effect of non-uniform fiber Imgth on bundle strength 
Number of fibers = 60 Number of trials = 500 



Mean Strength, 
Xy(MPa) 

cav 

«* 

P*(3£Pa) 

0.0 

630.1 

0.0438 

23.79 

643.3 

0.5 

615.8 

0.0442 

24.51 

628.7 

1.0 

573.8 

0.0464 

22.36 

586.5 

2.0 

471.4 

0.0689 

15.61 

486.4 

3.0 

385.0 

0.0850 

12.68 

399.9 



Figure 2.4: Effect of non-uniform fiber lengths on bundle strength. 

1 
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TABLE 2.4 


Modified bundle strength predictions for 


“/ 

Mean Strength (MPa) 

Coefficient of Variation 

8.2 

658. 

0.045 

6.0 

605. 

0.053 

4.0 

533. 

0.067 

2.0 

416. 

0.100 

experiment: 




682. 

0.220 


increased moderately and the mean strength decreased for a decreasing value in the shape 
parameter. Complete results are given in Table 2.4. 


SUMMARY 


The modified bundle theory, even with the irregularities which have been introduced 
in the previous section, does not adequately predict the variation which is seen experimen- 
tally for the post-matrix-cracldng longitudinal strength of a layer, X^. This suggests that 
the model is ignoring one or more important mechanisms. One possible area to be focused 
on is the fiber/matrix interface; the interfacial shear stress is not truly a constant, stress 
concentrations are present on the fiber at the crack face, also the matrix environment 
makes the reloading of failed fibers away from the failure location possible. However, 
while there are problems with the modified bundle analysis in its present form, the results 
of the analyses which simply allow for a complete all-at-once failure in the longitudinal 
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direction are genraally good. This is encouraging for the basic modal analysis and load 
redistribution scheme. 
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Generalized 
Shear-Lag Analysis 


The onset of transverse matrix cracking has been found to be a key occurrence in the 
laminat e failure process, and as a result much research has been devoted to its modeling. 
Continuum damage models (Talreja, 1985; Nuismer and Tan, 1988) as well as many shear- 
lag approaches (Garret and Bailey, 1977; Reifsnid^, 1978; Laws and Dvorak, 1988; Lee 
and Danid, 1991) have been proposed for analyzing transverse cracking in cross-ply 
laminate s, but very little work has addressed the problem for layers of arbitrary 
orientation. In this chapter, the method developed by Lee and Daniel for determining the 
shear-lag parameters is reviewed and extended. Their original ^plication of the method 
was for a cross-ply [®m/^n], lan^de; here, the shear-lag relationship will be developed 
for a generalized system comprised of three arbitrarily oriented layers of the form 
[e/4>/t|r],. This generalized shear-lag analysis will then be used in Chapter 4 to aid in the 
solution of the elasticity problem for the region of a general symmetric laminate between 
two parallel transverse matrix cracks. 
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Figure 3.1: Traditional shear-lag model. 
SHEAR-LAG METHODOLOGY 


Shear-lag analyses are generally set up to relate the shear stress present at the interface 
between two contiguous bodies to the differaice in their respective average displacements. 
That is to say, for the bodies shown in Figure 3.1, a shear-lag rule applied in the 
traditional manner would define the shear stress presoit at the interface between 
bodies 1 and 2 to be proportional to the difference in the av^ge displacements ( and 
^ observed in each of these bodies. 

x^=H[u^-u^) ( 3 . 1 ) 

The proportionality constant Ff in the above equation is referred to as the shear-lag 
parameter, and is a property of the entire system. 

In regard to the application of shear-lag methods to composite laminates, layers 1 and 
2 of Figure 3.1 generally refer to two orthogonal layers which comprise one half of a 
symmetric laminate system. One of these layers is considered to have undergone some 


22 



Chapter 3* CjcncrdlizGd An&lysis 


form of damage, typically transverse matrix cracking, which results in the displacement 
variation between the layers. If it were not for the damage being present, laminated piatp 
theory would predict a uniform displacement field through the thickness of the system! - 
in the absence of bending moments and non-symmetric configurations. Various methods 
exist to evaluate H, employing mechanics, energy principles, hindsight, and combinations 
of the above. The most methodical and robust approach, in this author’s opinion, was 
proposed by Lee and Daniel (1991) and expanded to two dimensions by Tsai, Daniel and 
Lee (1990). This is the method which will be adopted for this report, and the remainder 
of this chapter is devoted to extending this method to a three layer system. 


MODEL DERIVATION 

The derivation is begun by expficidy defining the system of interest. The laminate is 
assumed to be symmetric, and therefore attention can be focused on only those layers lying 
on one side of the laminate’s mid-surface. Thus, a six layer laminate reduces to only a 
three layer system. Considering the case where the middle layer of this three layer system 
has experienced transverse cracking, the overaU goal for this chapter and the next is to 
determine the effects of this damage on the layer’s behavior. That is to say, that the 
nuddle layer is the main focus of interest, and therefore, the definition of the problem is 
setup with this layer in mind. The layer numbering scheme is established in a manner such 
that the damaged layer is referred to as layer 1, and the two adjacent layers as layers 2 and 
3. This is shown in Figure 3.2. 

The solution of the problem will be conducted in the material coordinate system of layer 
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1. This requires a transformation from the global x -y coordinates of the laminate to the 
x'-y' coordinates of layer 1 (note that this is an orthogonal transformation about the z 
axis such that z'= z). The system is now suffidoitly defined for the derivation to proceed. 

Key to this derivation is the assumption that the out-of-plane shear stresses vary linearly 
through the thickness within each layer. As a result, a quadratic displacement field 
through the thickness direction for each of the layers is proposed 


u\{x\y\z)= 

(3.2) 

v\{x\y\z) = ^4 z* + fljZ + «6 

(3.3) 

u = a, z^ + a,z + 0, 

(3.4) 

v\{x\y\£i = + "i2 

(3.5) 

«'3(x',y',z) = +«i4Z + «i5 

(3.6) 

v'3(x;y',z) = ai6Z^+«i7Z+"i« 

(3.7) 


with u\ and v', designating displacement in the x’ and y’ directions, respectively, for 
layer i (where i = 1,2,3), and the coefficients representing undetermined 
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functions of x' and y'. 

The constitutive relation for transverse stresses and strains for layer i is given below. 


fM _ 

Q'ss 

Q'as 



1 

P'as 


i 

[v.. 


(3.8) 


^o summation is implied by the repeated subscript in equation (3.8)). The primes on the 
stiffness terms are present to signify that they have been transformed to the x'-y' 
coordinate system. These equations can be considm^ed separately from the in-plane 
constitutive equations because of the assumed orthotropy of the layers (Reismann/Wether- 
hold, 1988). 

The appropriate strain-displacement relations are given in equation (3.9). 


yO) _ ^ 

dz 


Y« = 

dz 


(3.9a,b) 


These equations neglect the influence of deformation in the out-of-plane direction, i.e. , 
and . Utilizing the displacement equations (3.2 - 3.7) as well as the strain-displacement 

dy' 

relations (3.9a, b) in the constitutive law yields the following system of equations relating 
the out-of-plane shear stresses to the unknown coefficients of the displacement equations. 

f2fliz+a2l 


f \ 






2a^ +flj 



' = [QV 
2 

20jZ+a^ 

Wz\ 

2aioZ+flii 


^ = [QV 

3 


k'J 

2a^^z*a„ 


(3.10) 


(3.11) 


(3.12) 
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Note that use has been made of the fact that the transformed coordinates are aligned with 


the material coordinates of layer 1 , and thus = [Q]j . 

Considering the boundary conditions (see Figure 3.3), states of zero shear exist on the 
top face of the laminate (z = h^ +A 2 +A 3 ) and at the mid-plane (z = 0 ), 




Ty^(Z=A,+*2+*3), 


2 



(3.13) 


V.(J=0)1 

i.H 


i J 


and conditions of continuity exist at the two layCT interfaces 


(3.14) 



+ 

II 


^0-2) 

Vz 


Ty.,(z=Ai+h 3 ) 

^ K 
1 

^0-2) 
[Vz J 


(3,15) 
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'*3)1 _ jV.(i =*3)l _ k”."| ^3 

M'-*3)J, lv.(«'*3)|, 

We have introduced the variable definition to designate the 

shear stresses present at the interfaces. These stresses are undetermined functions of x' 
and y'. Applying the above boundary conditions to equations (3.10 - 3.12), the 
coefficients for the linear and quadratic z terms can be solved for as functions of the 
interfadal shear stresses, and layer stiffnesses and thicknesses. This results in displace- 
ment equations containing the shear stresses and the "a" coefficients of the zero-order z 
terms as unknowns. 




Requiring displacement continuity for u' and v' at the (1-2) and (1-3) interfaces yields 
expressions for the last remaining unknown coefficients, 
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where 


= [QTxRA 




M-2) 

.( 1 - 2 ) 

‘y's 


t + 





(3.21) 


ih,^iH){hrh,) 

2h, 


Ra 


Ji. 

2h, 




{hi*h,){h,*2h^*h,) 

2*2 




fh-^h,*h,)h, 

2h, 



Equations (3.20) and (3.21) represent a system of four independent equations in terms of 
six unknowns (flj, a^, a^, dnd nj,). There is not sufficient information available 

to ejqilidtly solve for each of the unknown terms. However, by averaging the displace- 
ment equations and subtracting (in order to find the difference in average displacements 
required for a shear lag analysis), all the unknown coefficients will drop out of the 
analysis. 

The average displacement within each of the layers can be determined by integrating 
the corresponding displacement equation (i.e., equations (3.17) through (3.19)) over z and 
dividing by the req>ective layer thickness. This procedure results in the following set of 


equations for the average through-the-thickness displacements (note: these displacements 


are still functions of x' and y'). 

fiT' 






i: 


( 3 . 22 ) 
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where 


u 


= [Q']2\ 


Jl-2) 

0 - 2 ) 


+ < 


*12 


fr:' 1 


0-3)1 

f 

^ 3 ] 

• = ■ 

Vz 

«15 


0-3)[^^® * 
^y'z ] 

«18 


(3.23) 


(3.24) 


«. = 





2h\ 

j[{*i ~(K ] -{*1 



R 


10 


^3 

6 


Subtracting to find the appropriate "differences" in the average di^lacements results in 


equations which are functions of the layer stiffness and thickness terms (which are known). 


f3-"i5 

K "«18 


"3-«9l 
re “^12 1 


and 


the interfacial shear stresses (which are unknown), and the expressions 

involving the remaining coefficients from the displacement equations (these are 
known from equations (3.20) and (3.21)). The final forms of these difference in 
displacement equations are given below: 



m 


^0-2) 

0 - 2 ) 

^y'z 


*[B]\ 



(3.25) 
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,0-2) 


_0-3) 

U tt 3 


Vz 

■»[£)]• 

Vz 

v'l-v'a] 

. = [cy 

,0-^ 

,0-3) 

yrz J 


(3.26) 


where 


[^1 - - y[<?']? 

[B] ' 

[Cl- ^[<?i;‘ 

[O] - y[<?i;‘ * y[<?']? 

At this point, as a check, the three layer model can be compared to Tsai, Daniel and 
Lee’s two layer model. If layer 3 were no longer preset, the (1-3) intefece would 
become the laminate’s mid-surface, and as a result the shear stresses andty^j^^ 

would be equal to zero due to symmetry considerations. Under these conditions, equation 
(3.25) becomes 


“r«2 

f 


II 

1 

v\-v\ 

1 


lOl'y * [O'E'y 



r 

,0-2) 

i 

Vz 


,0-2) 


[v* J 


which agrees with the published solution (Tsai, et al., 1990) for a two layer system. 


Equations (3.25) and (3.26) can be inv«ted to solve for the shear stresses as a function 
of the difference in the average in-plane dii^lacements for each of the layers 


,0-2) 

Vz 


** 1 ® 2 

■ .[J]. 

« l“ K 3 

,0-2) 

J 

• • [»]■ 


1- - 


(3.26) 
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^(1-3) 


r ' 


r— p- — ^ 

T , 

■ -W' 

U 1 tf 2 


“l-»3 

[S'z J 


■Mi]' 



where 


(3.27) 


[t] = [[D]-(C] [.<]-■[«]]■* (3-28) 

[ff] - [.4]-‘*[.4]-‘[B][L][C][-4]- (3.29) 

[7] - (3.30) 

m - -[i][C)(^]-‘ (3.31) 


In this form, resemblance to traditional shear-lag analyses can be seen, with the out-of- 
plane shear stresses at the interfaces being stated explicitly in terms of the differences in 
the average in-plane displacements. The matrices [H], [J], [J^ and [L] are the shear-lag 

parameters for the system. These parameters are constants which are completely defined 
by the equations above, and contain only elastic constants and lay^ thicknesses. 


SUMMARY 


Li this chapter a generalized shear-lag analysis has been developed to model a system 
comprised of three arbitrarily oriented layers (and thus, a six layer symmetric laminate). 
The only assumptions made in the development of this model were the linear variation in 
the out-of-plane shear stresses through the thickness of each layer, and the form of the 
stress-strain constitutive relationship. These assumptions along with stress and displace- 
ment continuity considerations then led to the shear-lag relationship of equations (3.26) and 
(3.27). This relationship shows the system to be coupled: the shear stresses at the interface 
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between layer 1 and layer 2 are a function of not only of the displacements in layers 1 and 
2, but also the displacement of layer 3. Similarly, the shear stress at the (1-3) interface 
is cou|>led to the displacements in all three layCTS. Were it not for a methodical approach 
to the derivation of the shear-lag relationship, this coupling effect might very easily be 
overlooked, and the functional form of the effect would most certainly not be otherwise 
intuitively determined. 
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Elasticity Solution 

of the 

Three-Layer Problem 


As discussed in the Introduction, the failure process in composite laminates often 
initiates with the onset of transverse matrix cracking within the individual layers. The 
initiation of transverse cracking within a layer occurs when the transverse stress present 
within the layer reaches the local strength. This strength value is typically referred to as 
the first matrix cracking strength. When dealing with brittle materials (e.g., ceramics) 
the purpose of designing with composite materials as opposed to monolithic materials is 
to reduce the flaw sensitivity of the material system. While this desired effect is realized 
in a rdative sense, brittle matrix composite materials still possess an inherent flaw 
sensitivity, and as a result the observed strength characteristics for the composite contain 
scatter. In this report the material strengths are considered to be random variables 
represented by WeibuU distributions. However, even after transverse cracking has 
occurred within individual layers, the system as a whole can continue to survive due to the 
in-situ constraints of the laminate system on these individual layers. In fact, due to 
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shearing effects between neighboring undamaged and damaged layers, the damaged layers 
themselves can still continue to carry load! In this chapter, the reloading of a damaged 
layer in the region between two parallel transverse matrix cracks is investigated in detail. 
This analysis will be carried out for a system of three arbitrarily oriented layers under 
general in-plane applied load conditions. It will be shown through homogenization 
techniques that this solution is sufficient to model the damage effects of transverse matrix 
cracks in any symmetric laminate. 

DEFINING THE PROBLEM 

The investigation begins by considering a goieral symmetric multilayered laminate 
having a total of 2n layers. The corresponding layer thicknesses and orientations are 
by hj and 0j, respectively, for / = 1 to 2n. It is assumed that the laminate is 
subjected to a gaieral in-plane load of (2P^f2Py,2P^) . Owing to the laminate’s 
symmetry, attention is restricted to only those layers above the mid-surface (i.e., layCTS 
1 through »), thereby reducing the complexity of the problem by a factor of two. As a 
result of the applied load, transverse matrix cracking begins to develop in one of the mnCT 
layers, namely, layer i . Figure 4. 1 depicts a schematic representation of the system just 

described. 

The goal is to determine how the formation of these transverse cracks effects the load- 
carrying capabilities of the damaged layer, and thus in turn the global level behavior of the 
entire laminate. To achieve this, the stress state for the region between two paraUel 
transverse cracks of the damaged layer needs to be determined at the local level. This 
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A 

z 



X 

F^re 4.1: Laminate undergoing transverse matrix craddng in layer i. 


stress state is a function of position in all three spatial coordinates, i.e. { o }, = { o (x,y,z)}, . 
However, because interested is focused on the in-plane effects of this damage (i.e., the 
problem deals with a thin plate within the analytical confines of classical laminatfvt plate 
theory), only functional relationships with respect to x and y are of main interest. Thus, 
the solution for the average through-the-thickness stress within the layer, {a}. = {a(x,y)}., 
is pursued. As a note, throughout the remainder of this chapter, reference to a variable 
as an average value or the use of the over bar symbol will always refer to a through-the- 
thickness average unless otherwise specifically stated. 

At the local perspective from which this problem is being approached, only the stress 
solution for the area in the damaged layer between the two transverse cracks is of concern. 
The knowledge of the stress state throughout the remainder of the laminate is only of 
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interest as it affects the damaged region. Thus, in an effort to reduce the complexity of 
the problem, it is proposed that all layers, other than the damaged layer of interest, be 
homogenized. This action results in a system which is elastically equivalent at the global 
level, and has been effectively reduced from n-layers to only 3-layers. This new 
homogenized system is shown in Figure 4.2. In the new system, the layer numbering 
scheme has been reordered such that layer 1 refers to the damaged layer and layers 2 and 3 
the two adjacent layers, with the tildes (-) used to differentiate between the 
homogenized and original systems. The elastic properties and layer thickness for the 
damaged layer in the new system, layer i , and in the original system, layer i , remain the 

same, i.e., 

[S]i = [S], (4-1) 

hi = h, (4.2) 

The properties for layers 2 and 3 are now giv«i by a thickness weighted avwage 
eqnessed as 

(4.3) 

E*/ 

Here /ij and fcj are the layer thicknesses in the new system, and [s ]j and [s ]§ are the 
compliance matrices for the homogenized layers. The [s ] notation signifies that the 
compliance matrices have been transformed from the local material coordinates of the layer 
to the global coordinates of the laminate. Since layers 2 and 3 possess the same elastic 
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Figure 4.2: Elastically equivalent three layer system. 


pnqierties in the homogenized system, the problem to be solved is further simplified. 

Figure 4.3(a) shows a top view of the laminate. The matrix cracks in layer 1 are 
designated by the dashed lines. It becomes apparent from this view of the problem that 
the natural coordinate system in which to pursue the solution is the local material 
coordinate system of the damaged layer. Thus, the problem is set up in the x'-y' 
coordinates of layer I . Note that this is an orthogonal transformation about the z axis 
such that z'= z. A detailed view of the representative system for this problem is shown 
in Figure 4.3(b). From this point on, all analysis will be conducted for the homogenized 
system. Therefore, no ambiguity should result from dropping the tildes ( - ) from the layer 
designations in order to ease notation. 
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Figure 4.3 (a) Top view of the laminate, (b) Detailed view of the representative 

volume element to be solved. 

PROBLEM SOLUTION 

The solution begins with a formulation of the in-plane equilibrium equations for the 
three-layer system. Ck)nsider blocks of diffo^ntial size in the in-plane directions and with 
finite thi ck nesse s corresponding to the individual layer thicknesses in the out-of-plane 
direction (see Figure 4.4). Summation of forces in the x' and y' directions for layer 1 
yields. 
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Figure 4.4; Free body diagrams for individual layers in a three-layer system. 
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52 fW ^ 0 = 
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(4.5a) 
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dx’ + - x'^^l^^dx'dy 
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a-3) a-2, . 
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(4.5b) 


= 0 = -N^WAN^^*^dy' 


a) 


dK 


a) 


ay' 


dx’-N^W 


xy- 


cbc^ 


dy' + 'dy 


(4.6a) 


0-» (l-» . 

■ dx’ dy’ 


(4.6b) 


In a <^miiar manner, the forces acting on layCTS 2 and 3 may be summed. The resulting 
local equilibrium equations for all of the layers are grouped together and presented in 
equations (4.7) through (4.12) below. 



,0-2) _ 

dN^ dN^i’ 

(4.7) 
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,0-2) _ 
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dN^^ dN^l^ 

(4.9) 
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^a-3) _ 



(4.11) 

X z 

at' 

ay' 

Jl-3) _ 


+ 

(4.12) 


ax' 

dy' 


We "lalrp. the assumption that the effects of the transverse cracking are direction-specific. 
That is to say, the influence of this form of damage is seen as the observer moves in a 
direction perpendicular to the matrix crack face (y' direction); and if the observer moves 
parallel to the fiber (x' direction), no damage influence is detected. Thus, all parameters 
are constant with respect to x' ( i.e., -?7=0). This simplifies the local equilibrium 
equations, resulting in ordinary differential equations with independent variable x' only. 


= 

dN^l- 

(4.13) 

Vz 

dy' 




(4.14) 
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dy' 


^(1-3) ^(1-2) _ 
Vz ~Vz " 

dNxY 

dy' 

(4.15) 

a-3)_ n-2) _ 

dN^^ 

dy' 

(4.16) 

.,0-3) _ 
“Vz " 

dy' 

(4.17) 

_ 1 - 0 - 3 ) _ 

dN^^ 

(4.18) 

^y'z 

dy' 



A further result of this assumption is that displacement continuity between each of the 
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layers in the x' direction is maintained. Thus, 

u\ = = tt 3 (4.19) 

The objective is to determine the stress state for layer 1. Therefore it is reasonable to 
look at equations (4.15) and (4.16), which rdate the interfacial shear stresses to the in- 
plane loads for layer 1. Begin with equation (4.16) which specifically deals with the rate 
change of normal loading in the y' direction as we move away from the crack face. 

- V. "V. 

This equation contains both in-plane effects (i.e., the loads) and out-of-plane effects (i.e., 
the shear stresses), but can be written entirely in terms of in-plane effects by making use 
of the shear-lag equations (3.26) and (3.27). Using these equations, the interfacial shear 
stresses appearing in this equilibrium equation can be written in terms of the in-plane 
displacements. 


,( 1 - 2 ) 

h'Z 


,( 1 - 3 ) 


^12 (“ l““.2) ^22('' 1 ''2) *^12(“ 1~“ 3) 


(4.20) 

(4.21) 


From the displacement continuity assumption, the difference in the displacement in the x' 
direction between each layer must be zero, i.e. = {“^““ 3 ) = 0- The local 

equilibrium equation thereby becomes, 

-fl=) (v'l - (Ln ) (V'. (4.22) 

DiffCTentiating this equation with respect to y', and applying the strain-displacement 
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relation 


V = 


dv' 

dy' 


(4.23) 


yields a governing differential equation containing both unknown load and strain terms. 


y 


dy 


- S®) * [hi (4.24) 


LAYER THICKNESS EFFECT 

Thus far the governing equation has been modified iirom one with mixed in- plan ^ and 
out-of-plane effects to an equation solely in terms of in-plane terms. The next step 
involves a further moditication to this equation such that the unknowns T^resent the same 
physical phenomena (i.e., either strains or loads). By utilizing the constitutive relation for 
the individual layers, the strains appearing in equation (4.24) can be expressed in terms of 
the in-plane loads. 



(4.25) 


(4.26) 


(4.27) 


It is important to note that though the only stresses contained in these equations are in- 
plane, this does not imply that a plane stress state is being assumed. Because the 
individual layers are orthotropic, the in-plane strain and out-of-plane stress are uncoupled. 
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Therefore, equations (4.25) through (4.27) do not preclude the existence of out-of-plane 
shear stresses; this is indeed important due to the role they play in the reloading scheme. 
The on the compliance toms are present to signify that they have been transformed 
to the x'-y' coordinate system (note that no transformations for the compliance terms for 
layor 1 are necessary). 

At this point, two specific examples (one for a two-layo: case and another for a three- 
layer case, see Figure 4.5) and an inconsistency is detected. In the two-layer case, a layer 
containing transverse cracks is faced on one side by an undamaged layer which provides 
for reloading the damaged layer. The related three-layer case has been created by basically 
"flitting up" the undamaged layer from the two-layer case into two separate layers and 
placing one on each fece of the damaged layer. The governing differential equation for 
the two-layer case is given below. 




dy‘ 




'2-layer east 


(4.28) 


f 2-layer ease 


where ^ designates the shear-lag parameter for the two layer case (Lee and Daniel, 1991); 
the governing equation for the three-layer case is given in equation (4.24). 

The elastic properties are the same for the two cases, therefore the loads are related as 


follows. 


(4.29) 




(4.30) 


Calculation of the strains using the constitutive relations above leads to the incorrect 
prediction that the difference between the strains of layer 1 and layer 2 are the same for 
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Figure 4.5: Related two-layer and three-layer cases. 


both cases. In actuality, it is known that there is a thickness effect present such that as the 
thickness of the facing layer decreases, the difference in the strains also decreases. Thus, 
the relationship sought to predict is: 

(«i) 


In orda to achieve this prediction, the introduction of a correction term becomes 
necessary. 

In proposing the form of the thickness correction terms, terms are defined such that the 
rate of reloading is the same in both cases. Thus the following relationship is required to 
hold true. 




l k -layer cast \ dy'^ h -layer case 




.'2 


This condition may be restated using the shear-lag relationships. 


(4.32) 




(4.33) 
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wh^ thickness correction terms, 62 ^3 > lay^ 2 and 3 respectively, have been 

introduced to modify the differrace in strain terms. For this example case, upon 
evaluating the shear-lag parameters, we find 

= 2xH (4.34) 

Additionally, because the elastic properties are the same in the two facing layers, the 
strains in these layers are the same. 



(4.35) 


Substituting this information into equation (4.33) results in the following condition on the 
thickness correction terms, 

8 , » *, = I (4.36) 


The ensuing additional conditions are imposed, 


62 = 63 = 7 

if 

II 

(4.37) 

62 = 0 ; 83 = 1 

if 

0 

II 

(4.38) 

«2 = i = 0 

if 

0 

II 

ptT 

(4.39) 


The following forms which satisfy the above conditions are proposed for the correction 


terms. 


6 - ^ . 5 - *3 

^ " 2(*2+A3) 


(4.40) 


An alternative interpretation (Tsai, 1993) to the discrepancy between the two-laya: and 
three-layer formulations is that a facing layer, regardless of its thickness (or lack thereof). 
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provides reloading support to the system; and thus the fact that the solution for predicted 
shear force transmitted to the damaged layer is discontinuous at the transition point from 
a two-layer to a three-layer system should not be unexpected. 

This author believes that "real life" lies somewhere in between these two interpreta- 
tions. That is to say, that a three-layer configuration provides a more efficient reloading 
rnf^rhanigm , yet there Still must be some thickness effect present. The perspective of Tsai 
would over-predict the reloading capability of the system, and thus err unconservatively. 
Whereas the interpretation of the thickness effect presented in this section approaches the 
correct solution in the limit, while erring on the conservative side. Thus, it is this 
thickness effect model which is adopted for the remainder of this body of work. 


SOLUTION (continued) 


In the homogenized system, the elastic parameters for layers 2 and 3 are the same, and 
as a result The governing differential equation for the three-layer case thus 

becomes 


dy'^ 



(4.41) 


where is a known constant which is a function of the shear-lag parameters and the 
layer thicknesses. 

Hn ' -^) " -^)1 (‘‘•“ 2 ) 

Substitution for the strains using the constitutive relations (equations 4.25-4.27) produces 
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a governing diffmoitial equation entirdy expressed in terms of in-plane loads. 



However, as it stands this one equation contains six unknowns and thus five additional 
independent equations must be found. From global equilibrium considerations (see 
Appendix A) three of the necessary equations are obtained. 



(4.44) 


(4.45) 

k- 

(4.46) 


These can be used to solve for the in-plane loads of layer 2 in terms of the in-plane loads 
of layer 1 and the loads Py and P^y (transformed from the known applied loads 
PjfiPy and P^, see Appendix B). Continuing the previously made assumption of 
diic pIareTnent continuity in the x' direction implies that the normal strains in this direction 
for each layer are equal. 



(4.47) 


This equation can be used to solve for in terms of the remaining unknown terms. 
Substitution of these results into equation (4.43) produces a second order, ordinary 
differential equation with only two remaining unknowns, and N^y. 


dy'^ 




(4.48) 
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The coefficients a , r\ and P are constants which are algebraic functions of the elastic and 
shear-lag parameters, and P also contains the applied loading terms, 


a 


( pW c'® ) 
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-1 
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3i2 •> 12 
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(4.49) 
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(4.50) 




(4.51) 


where 


A. = 
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1 

A2+A3 

1 

A2+A3 

1 


^ 12 11 

O 22 12 


«p(2) 

^ 26 16 


' cW VV S'^'^ 'i 

h^ 

'si? si?r[»g. ' 

^ *, **2<-*3) i *l *2^*3, 

cjO) C'O) W 


(4.52a) 


(4.52b) 


(4.52c) 


The last required equation comes from the set of local equilibrium equations. Using 
equation (4.15), and applying the shear-lag relations along with the necessary thickness 
correction terms, a second order ordinary differential equation for I^y is achieved which 
is similar in form to equation (4.41). 


dy'^ 



(4.53) 


49 


Chapter 4: Elasticity Solutioa 


where 


* 2(^2 + Aj) ~^)] (4.54) 


Following a procedure similar to the one outlined above, this differential equation can be 
reformulated such that it contains the known constants ilrjj , a , t) and p , and the 
unknown loads and N^y. 




(4.55) 


Thus with equations (4.48) and (4.55) a system of two coupled, second order, ordinary 
differmtial equations is obtained. The boundary conditions for this problem are homoge- 
neous owing to the traction-free crack surfaces, 

W,?’(y'=0) = Ar“(y'=2*) = 0 (4.56) 

= wi?.(j>'-2») - 0 (4.57) 

where 2s is the distance between two paralld matrix cracks. This system can be 

uncoupled at the e^)ense of increasing the order of differentiation, resulting in the 

following fourth order, homogeneous differential equation, 

, d%<!> 

=0 (4.58) 

with the new boundary conditions, 

A^?^0) = N^\2s) = 0 (4.59) 
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and 


d^N^}\0) d^N^^(2s) 

dy'^ dy'^ 


P ^^22 


(4.60) 


The solution of this equation for Ny^{y') then leads directly to N^^{y') and N^y{y') 
being known as well. The final forms of the solution for all three in-plane loads of layer 
1 can now be written as 


Ny\y') = BjSinh(Ay') +B2CO^{ky') +B^ 

(4.61) 

Ny\y') = B^sinh(A,y') +H 5 Codi(Xy')+B 5 

(4.62) 

Nyyiy') = .B7sinh(Ay') +BgCosh(>.y')+^9 

(4.63) 


with e7q>ressions for k and the coefficients through given below. 


k = +y^at|f22 + niki2 
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This solution for cases where a full goieral in-plane loading is applied is dependent upon 
a non-zero value for the shear-lag terms ijr-jj and i|r| 2 . The value of ijfjj will always be 
non-zero; however, for certain laminate configurations i|rj 2 = 0. Further details are 
discussed in Appendix C. 

F inally , the shear stresses at the interfaces can be calculated. From the local 
equilibrium equations, the shear stresses at the (1-2) interface are known in terms of the 
in-plane loads of layer 2, 


,a-2) _ 


•xz 


ds’ 



dy' 


(4.64) 

(4.65) 


Combining the solutions for the in-plane loads of layer 1 with the overall equilibrium 
equations (4.44) through (4.46), the loads for layer 2 may be determined and equations 
(4.64) and (4.65) can be evaluated. 

^a-2) ^ _!!!L-^-RjXco^{Xy')-B^Xsarii{Xy')] (4.66) 

= -^[-ll4Xcosh(Xy')-BjXsinh(Xy')] (4.67) 

H2+/I3 
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Similarly, from the local equilibrium equations for layer 1, 


_a-3) 

__0-2) _ 

dN% 

(4.68) 



dy' 


^a-3) 

S'* 

-d-2) _ 

^y'x 

dy' 

(4.69) 


These equations yield expressions for the shear stresses at the (1-3) interface. 

=[i ^]rB7Xcosh(Xy')+fl,A.sinh(Xy')] (4.70) 

= [l--^][B^Xcosh(A.y0+^5^“*(^y'^^ (4.71) 

EFFECTIVE ELASTIC BEHAVIOR 

The in-plane loads of a damaged layer within the region between two transverse matrix 
cracks have now been completely solved for. The next step is to determine how this layer 
will effectively behave, from a macro-level viewpoint, within the laminate. In examining 
this bdiavior, the interest is in the in-plane-average req>onses (in contrast to the average 
through-the-thickness that has been discussed up to this point). In order to emphasize that 
the averages being considered in this section are with regard to both the thickness of the 
layer as well as the in-plane length, a double over bar notation (*) is introduced to signify 
that the average is with respect to both dimensions, and the single over bar notation is 
continued to denote averages taken with respect to thickness only. 

The effective compliance of layer 1, [^]i, is defined by the relationship between the 
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in-plane-average strain of the layer» \€/| , and the in-plane-average stress, \a/| , as shown 
below, 


II 


^11 ® 
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^12 4 0 
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1 

. Vjr 'j 


(4.72) 


Modding the material in this fashion provides a measure of the effective "secant” (as 
opposed to "tangential”) behavior of the lay». Ihat is, if the applied strain is assumed to 
be removed from the layer, the layer’s load will decrease to zero (the origin of the stress- 
strain diagram) linearly. Upon reloading, the lay^ will again behave linearly until 
reaching the previous load state. Further straining beyond this point can result in 
additional damage accumulation and a reduction in the effective modulus as just described. 
Expressions for the in-plane-average stress can be determined by integrating the stress 
equations (found by dividing equations 4.61 through 4.63 for the in-plane loads by the 
layer thickness; i.e., over the crack bounded r^on. Note that these 

equations are independent of x'; as a result, averaging is only necessary over y'. For the 
average stress of layer 1 in the x' direction, this yields 


=0) 


2s J 


(4.73) 






[cosh(2A5)-l] 


^2 ^3 

sinh(2A.s) + — 


and similarly for the remaining stress terms. 
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= — -i-[cosh(2Xs)-l] + -i-siiih(2Xs) + ^ (4.74) 

^ 2k^Xs A, 

f], = _iL_[cosh(2Xs)-l] . -^siiih(2Xrt . ^ (‘‘■75) 

" 2*,Xj^ 2»,Xs *, 

The above equations can be used to determine the average in-plane stresses based on the 
fer-field applied loads (which the coefficients through depraid on) and the crack 
spacing (2s) . The size of this crack spacing can be determined in an iterative manner. 
Transverse cracking occurs when the in-plane transverse stress reaches the material s 
transverse strength (Xj). From the form of equation (4.62), it can be seen that the 
transverse stress reaches its maximum value at the midpoint of the crack spacing. 
Therefore, by iterating on the value of the half crack spacing, s , until the calculated value 
of the stress at the midpoint reaches the strength value, i.e. o^‘iy'=s) = the crack 
spacing can be determined. 

Remembering that it is the effective global behavior of the layer within the laminate that 
is being modeled, it is reasonable to assume that the average strains appearing in equation 
(4.72) are equal to the mid-plane strains of the laminate. 

(4.76) 

With this assumption, equation (4.72) relates the in-plane average stress of the layer (which 
can be calculated) to the strains of the laminate (which can either be measured or 
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calculated), all of which are known. The only unknowns are the four compliance terms. 

The assumption has previously been made that the damage effects due to transverse 
matrix cracking are restricted to the y' direction. This implies that the effective Young’s 
modulus in the x' direction should remain unchanged, 

- S, (4.77) 

and as a lesult the effisctive compliance tenn can be solved for. 

(4.78) 

This leaves only three remaining unknown terms in die effective compliance matrix, and 
these can be found using the three equations comprising equation (4.72). 



Finally, the effective engineering constants (i.e.. Young’s moduli, shear moduli, 
Poisson’s ratios - all being described in the secant manner) for the damaged layer can be 
determined from the effective compliance matrix. The definitions for the compliance terms 
as functions of the engineering constants are given below. 
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Substituting the effective compliance terms into the above relations, and solving for the 
aigineering constants yields 









These engineering constants describe how the damaged layer effectively behaves within the 
laminfltft system. In the next section, the ability to mod^ the effective behavior of the 
damaged material will be used to perform progressive fEiilure studies for example 
laminates. 


EXAMPLES 


In this section, specific example problems are investigated to verify that the generalized 
shear-lag model and the traditional shear-lag model agree for simple cross-ply systems, as 
well as to look at the more complex material and load configurations which the generalized 
model is capable of examining. The compute: code (PFRAC, see Appendix D) used in 
the analyses models the progressive failure of composite laminates under strain-controlled 
conditions. The material strengths are considered to be random variables of known distri- 
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TABLE 4.1 

Material Properties for Graphite/Epoxy 
(Lee and Daniel, 1990) 


Langitudinal modulus 

144.8 GPa 

Transverse modulus 

10.7 GPa 

In-plane shear modulus 

7.2 GPa 

Out-of-plane shear modulus 

3.8 GPa 

iu-plane major Poisson's ratio 

0.285 

Longitudinal tensile strength 

2167.7 MPa 

Transverse tensile strength 

54.5 MPa 

Longitudinal con^nessive strength 

-1440.3 MPa 

Transverse conq>i:es8ive strength 

-227.5 MPa 

In-plane shear strength 

82.0 MPa 

Interlaminar shear strength* 

113.0 MPa 

Ply thickness 

0.000127 m 


*Agarwil and Broutmui, 1980 


bution, and multiple modes of failure (i.e., longitudinal, transverse, and shear) are 
considered within each layer. As individual modes &il, correqx>nding stiffness reductions 
are made via the methods of the previous section as wdl as by traditional ply drop-off (or 
in this case mode drop-off) techniques. Monte Carlo methods are employed to determine 
the cumulative distribution function for the laminate. 

In the first example, agreement between a conventional shear-lag model (Lee and 
Daniel, 1991; Tsai, et al., 1990) and the gmeralized modd devdoped in Chapters 3 and 
4 is demonstrated for simple cross-ply laminates. The above referenced work examined 
a [0/90^1 graphite/epoxy (Gr/Ep) laminate. The material strengths were assumed to be 
deterministic, and thus the four 90° layers were considered to behave as one, resulting in 
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the conventional two-layer system. The material parameters can be foimd in Table 4.1. 
The Lee and Daniel analysis was performed for a load-controlled system, as opposed to 
the strain-controlled system mentioned above for the PFRAC code. Therefore, a modified 
version of the code is used for this example problem to allow for a direct comparison. 

This material system is examined for two laminate configurations, the [ 0 /^ 4 ]^ 
laminate studied by Lee and Daniel, and additionally, a laminate which is 

»q)lored to review the three layer capability of the generalized model. The analysis 
simulates a monotonic loading in the 0® direction. The [ 0 /^ 4 ]^ laminate was analyzed 
using both the generalized and the conventional shear-lag models, and the[0^/90^/0^J^ 
laminate was analyzed using the generalized model. The resulting stiess/strain curves are 
shown in Figure 4.6. As can be seen by the nearly coincident curves, very good 
agreement between the models is found. Additionally, it can be noted that a higher load 
to failure is predicted for the laminate. This can be attributed to the three- 

layer system’s generating lower interlaminar shear stresses than the two-layer system, and 
thus allowing efficient load tmsfer up to a higher applied load level. 

For the re mainin g examples, a four problem test matrix is established by examining two 
different material geometries as well as two different material systems (Gr/Ep and 
SiC/RBSN). The ability of the generalized shear-lag model to handle more complex 
laminate conEgurations is explored using laminate lay-ups of [0/15/80/-15/-80]^ and 
[0/30/60/-30/-60],. In order to perform the reliability analysis for the Gr/Ep material, 
Weibull strength parameters are needed, and are listed in Table 4.2 (note that it is assumed 
that the tensile and compressive strengths are equivalent). Experimentally determined 
values for the out-of-plane material properties of the SiC/RBSN system were not available. 
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Cjc 

Figure 4.6: Predicted stress-strain behavior for [0/904]* and [0w/904/0i/2]s laminates. 


For demonstration purposes, values of 130 MPa for the interlaminar shear strength and 
27.5 GPa for the out-of-plane shear modulus were assumed. 

The example problems will be conducted for cases of applied longitudinal strain, 
[€,] with free-edge conditions in the y direction (see Appoidix D). Even for this 
relatively simple applied strain state, the off-axis layers will still be e^riencing states of 
full in-plane loading. Therefore, examining the occurrences of transverse matrix cracking 
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( Ex y applied 

Figure 4.7: Realization of the stress-strain behavior for a [0/15/80/-15/-80]s 
Gr/Ep laminate. 


TABLE 4.2 

Weibull Strength Parameters for Gr/Ep 
(Wetherhold, 1986) 


a 


Longitudinal strength 

25 

1516.8 MPa 

Transverse strength 

10 

51.7 MPa 

In-plane ^ear strength 

15 

68.9 MPa 
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Figure 4.8: Realization of the stress-strain behavior for a [0/30/60/-30/-60]i 
Gr/Ep laminate. 

in any of the off-axis layers will require taking full advantage of the generalized shear-lag 
model. 

In the first set of figures, Figures 4.7 through 4.10, characteristic stress-strain 
re^nses are presented for each of the example problems. Each of these response curves 
are for one set of realizations of the material strengths. Contained within each figure are 
two curves comparing the results of the shear-lag model and the mode drop-off model. 
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Figure 4.9: Realization of the stress-strain behavior for a [0/15/80/-1 5/-80]s 
SiC/RBSN laminate. 

The algorithms for generating each of the curves use a drop-off technique for modeling 
stifhiess reductions due to longitudinal and shear ffdlures. However, for transverse failures 
(i.e., transverse matrix cracking), one algorithm uses the shear-lag model to determine the 
effective elastic properties of the damaged layw, and the other algorithm uses the drop-off 
method and removes that mode from future computations. In gaieral, both the shear-lag 
and mode drop-off techniques produced comparable results, eq>ecially for the Gr/Ep 
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Figure 4.10: Realization of the stress-strain behavior for a [0/30/60/-30/-60]s 
SiC/RBSN laminate. 

analyses. This may be partially due to the fact that, to avoid singularit y problems with the 
computer program, the stiffnesses were not exactly zero when using the drop-off technique. 
Thus, a damaged mode retained some load bearing capacity in both mcxlds. The stress- 
strain curves for all four examples demonstrate the evolution of damage in the material and 
die corresponding change in the global level cxinstitutive bdiavior. 

Figures 4.11 through 4. 14 contain characteristic plots of the crack density, p [ = (2s)"^ , 
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F^iire 4.11: Crack density as a function of ^jplied longitudinal strain for 
realization of [0/15/80/-15/-80]*Gi/Ep laminate. 

wh^ 2s is the crack spacing ] , in the damaged layers as a function of the applied strain. 
Note once again, these plots are for a particular realization of material strengths. For the 
Gr/£p mtamples, it is only the far off-axis layers ( ± 80 and ±60 ) in each case which 
experience transverse matrix cracking. In the SiC/RBSN laminates, cracking is seen in the 
15* and 30* near off-axis layers in addition to those further off-axis. The curves for all 
four examples are similar in nature, with increasing density and decreasing slope associated 
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( £v )^ppUed 

F^[iire 4.12: Crack density as a function of applied longitudinal strain for 
realization of [0/30/60/-30/-60]t Gi/Ep laminate. 

with increasing strain. Plateau r^ons ate present in some of the curves, and these can 
be attributed to times when there was a sufficient oiough decrease in the far-field stress, 
due to the failure of various internal modes, that the stress transferred into the damaged 
layers was not great enough to cause additional cracking. 

The last set of figures. Figures 4.15 through 4.18, contain the cumulative distribution 
curves for each of the example problems, with the probability of failure plotted as a 
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Figure 4.13: Crack density as a function of applied longitudinal strain for 
realization of [0/15/80/-15/-80]iSiC/RBSN laminate. 

function of the applied strain. The curves for the Gr/Ep laminates have the typical "S" 
shape common to cumulative distribution curves. The SiC/RBSN curves show a very large 
scatter in the ultimate strain (i.e., the strain to failure), with the maximum failure strain 
occurrence having a magnitude of approximately four times that of the minimum. The 
distribution curve for the failure of the [0/15/80/-15/-80], SiC/RBSN example in Figure 
4.17 levds off at approximately 0.003 strain and then b^ins increasing again at 0.0035 
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Figure 4.14: Crack density as a function of applied longitudinal strain for 
realization of [0/30/60/-30/-60]t SiC/RBSN lansinate. 

strain. This is characteristic of a bi-modal distribution. Furth«' examination of the Monte 
Carlo results showed that in fact this is so, with the lower aid of the distribution curve 
representative of those laminate realizations which failed due to the stiffness matrix 
violating the positive definiteness requirement, and the higher end of the distribution curve 
represoitative of those laminate realizations which failed due to a sudden 50% drop in the 
load carried by the laminate (these two failure critoia are aqilained in greater detail in 
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( Bx )i^pUed 

Figure 4.15: Cumulative distribution curve for [0/15/80/-15/-80]i Gr/Ep laminate. 


Appoidix D). 


SUMMARY 


A generalized shear-lag model has been derived to determine the average through-the- 
thickness stress state present in a layer undergoing transverse matrix cracking. The model 
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( Cx )i9pUed 

Figure 4.16: Cumulative distribution curve for [0/30/60/-30/-60]* Gr/Ep laminate. 


is capable of considering cracking in layers of arbitrary orientation, states of general in- 
plane applied loading, and laminates having a general symmetric stacking sequence. The 
model has been shown to agree with a conventional model for the case of a simple cross- 
ply laminate . Agreement was also found for the slightly more goieral three layer cross-ply 
system. Example problems were carried out for more graeral laminate configurations, 
namdy [0/15/80/-15/-80], and [0/30/60/-30/-60]^ laminates, with results showing the 
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Figure 4.17: Cumulative distribution curve for [0/15/80/-15/-80]s SiC/RBSN laimnate. 


evolution of damage in the material and the corresponding change in the global level 
constitutive behavior of the laminate. 
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Figure 4.18: Cumulative distribution curve for [0/30/60/-30/-60]* SiQRBSN laminate. 
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Location of Next 
Transverse Crack in 
Crack Bounded Region 


In the solution of Chapter 4, it was implicitly assumed that during damage progression, 
the transverse matrix cracks occurred at regularly ^ced intervals. This assumption 
allowed the solution of the stress state of the damaged layer to be reduced to the solution 
of a characteristic volume bounded by two transverse cracks and having a length equal to 
the crack spacing (2s). In this chapter the validity of this assumption will be investigated 
qualitatively by examining the probability density function for transverse crack location. 

In a classic paper by Oh and Finney (1970), the traditional Weibull analysis was 
extended to model the characteristics of the location of failure for brittle materials. It was 
demonstrated that for cases where the stress state present in the solid can be stated as a 
function of a far field reference stress as well as location, the probability density function 
^>df) for the location of failure can be determined analytically. Li this chapter, 
Wetherhold’s (1991) treatment of Oh and Finney’s work will be followed, and combined 
with the solution from Chapter 4 for the stress state present in the region between two 
transverse matrix cracks to arrive at the pdf for the next crack (i.e., failure) location. 
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STRESS SOLUTION 


Assuming that failure due to transverse cracking can be modeled using a probabilistic 
version of the maximum stress criterion (equation (2.3)), the critical stress which must be 
considered is the in-plane transverse stress. In Ch£q>ter 4, the solution was presented 
(equation (4.61-4.63)) for the state of in-plane stress of a damaged layer for the region 
between two transverse matrix cracks. As opposed to the goieral laminate configuration 
studied in Chapter 4, the laminate studied here is simplified for heuristic purposes. A 
uniaxially loaded, two layCT cross-ply configuration (see Figure 5.1.) is investigated. 
Corresponding to this simplified system, several reductions can be made to the general 
solution of Chapter 4. Due to the cross-ply configuration, the term becomes 
identically zero (see Appendix C). Because this is only a two layer problem, the thickness 
of the third layer is obviously zero, and drops out of all of tiie calculations. The uniaxial 
inatfing condition results in only the Pj, load (Py in the transformed coordinates) being 
presort. Taking the above effects into account and dividing the in-plane load (equation 
(4.62)) by the layer thickness, the average through-the-thickness transverse stress for the 
damaged layer can be written as 



A,Py 

1 -cosh(2Xs) 


sinh(2X^) 


sinh(Xy') + codi(Xy') - 1 


(5.1) 


The functional dependence on the location variable y' for this equation can be separated 
and defined by a new function r , 

This separation of variables allows the average stress to be erqnessed as a function of the 
applied load and the position. 
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CRACK LOCATION 


Consider that the region between the two transverse cracks consists of n elements. 
Each element may be uniquely specified by its position, y'^, and its volume, Ay'., where 
i = The probability that element i fails under the transvoiae stress present on the 

elemait is defined as (Oh and Finney, 1970; Wetherhold, 1991) 

/^[transverse failure of element i] s (5.4) 

From equation (5.3), the stress present on the element is a function of the position and the 
applied load, therefore 
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ft Ay i) (5.5) 

If the stress is considered to be constsnt over the volume of the element, then this equation 
can be restated using a failure density function. 

The function ijr, is a measure of the probability of feiilure per unit volume, and is a 
function of the reference load and the position. 

The conditional probability of failure for an element can also be studied. The 
probability that element Ay; fails on the load interval (P/, Py*dPy) given that the element 
has already survived a reference load of Py- can be expressed using the foUowing event 
statemoit, 

Pr[Ay'f fails for load 6 (Py',P/+dP/) | Ay; survivedP,*] = 

(5.7) 

Pr[(Ay; fails for load €{Py’,Py*-dPy)) f) (Ay; survivedPy>)] 

Pr[Ay; survived/^'] 

The probability that Ay; survived a load of /y- is a subset of the event that it fails on the 
interval [Py, Py- * dPy) . Thus, if the condition that the dement foils for a load e (Py-, Py + dPy) 
is met, the condition that it survived load Py' will automatically have been satisfied. 
Therefore, the conditional probability in equation (5.7) takes on the ensuing form. 

Pr[ Ay ; fails e (Py. Py+dPy) | Ay; survived Py'] 

Pr[Ay; fails c (J^'.iy'+dJ^')] 

Pr[Ay; survived Py^ 
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BG 


Ay; 


dPy 


*- dPy 


1 - 


t> ^yt) 


(5.8) 


The failures of each of the n elements comprising the entire damaged region are 
assumed to be independent events. Thus, the probability of survival for the whole region 
is givoi as the product of the individual survival probabilities. Combining this with 
equation (5.8) allows the detamination of the probability that the next failure in the region 
occurs at element AyJ for the interval (Py,Py+dPyy 


Pr[ Ay J fails H whole region survives Py'^ 


= dly' — 

aft- 1 


(5.9) 


n (> - Ay;)) 

b.y\dPy 


dP, 


1 - 0^,.{P,:y'„Ay-,) 


As the number of elements is allowed to increase, with n approaching infinity in the limit 
and the norm of the elemental volumes approaching zero, the finite volume AyJ can be 
replaced with the differential sized volume dy'.. Additionally, the denominator of equation 
(5.9) approaches one. 



and tile numwator, via the classical first order expansion employed by Weibull, can be 
written in the limit as having the following eiqionential form 
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Thus, the probability of equation (5.9) can be written as, 


Pr[demeiit dy' fails € Py+dPy'^ (5.10) 

n Tt^le region survives J^'] = hfPyty') dy' dPy 

where h is the joint density function for variables y' and Py', and is given below. 


h{Py,y') = 

Integrating equation (5.11), either with respect to Py or y' yields the marginal density 
functions for the location of failure, <l>(y'). or for the reference load at failure, g{Py), 
respectively. 

*(y-) ‘ f /‘{Pr.yyPy 

iV- 

g(P, } - f h(P,-,y'}dy' (5.13) 

y 

In evaluating the joint density function, it is assumed that the probability of failure per 
unit volume has a Weibull distribution. 




(5.14) 


where o. is the Weibull scale parameter and m is the Wdbull modulus. The following 
variable substitution is introduced in order to normalize the location variable. 
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y'e[0,2s] - 1 =^ ; 5e[0,l] (5.15) 

2s 

The joint density function then becomes, 




fPy') 

JB-1 


JH 

J 

r*(5)e^ 

-2sq\ — 



(5.16) 


where the constant q has been introduced to rq>resent the following integral, 

« . (5.17) 


The marginal density function for location which is being sought is then found by 
substituting expression (5.16) into equation (5.12), and integrating over the full range of 
the loading spectrum, i.e. (0,«»). 

= ( 5 . 18 ) 


By definition, if integration is carried out over all possible values of 5, thereby 
guaranteeing that all possible outcomes of <|) have been considered, this integral must be 
equal to one (this is a requirement of a pdf). Integration of equation (5.17) yields 

/>(U« - ^ - 1 (5.19) 

The density function must therefore be normalized in order to insure that the integral 
equals one. Denoting the normalized density function as the final form for the 
marginal doisity function for location of the next transverse matrix crack between two 
existing matrix cracks takes the following form 
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9 


(5.20) 


EXAMPLES 


In this section, the probability d^sity functions describing the location of the next 
transverse crack are examined for two mat^ial systems, the [ 02 / 962 ]^ SiC/RBSN laminate 
presented in Chapter 2 (Table 2.2) and die [0/90^]^ Gr/Ep laminate described in Chapter 
4 (Table 4. 1). For each of these cases, pdfs are gen^ated for a range of crack spacings. 
The evaluation of the constant q (equation (5.17)) contained in the final form of the 
solution requires numerical integration; all other calculations are capable of being 
determined in closed form. 

Figure 5.2 contains the results of the (normalized) pdf, ^ , plotted as a function of the 
location, ( , for the SiC/RBSN system. Five pdfs are dqiicted, each corresponding to a 
selected value of the crack spacing, 2s . These values range from a minimum of 2.5 mm 
to a maximum of 50 mm. Results for the Gr/Ep system, analogous to those just described 
for SiC/RBSN, are presented in Figure 5.3. 

Similar results are found for both material systems. The pdfs calculated for large 
crack spacings show large central regions of equal probability of failure. This corresponds 
to these regions being far enough away from the crack face to have uniformly reloaded. 
However, as the crack spacing gets smaller, the distribution narrows. This demonstrates 
the tendency towards cracking in the more highly stressed central regions. Consequently, 
in the early stages of damage, when the matrix crack density is relatively low, it appears 
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Figure 5 . 2 : pdf's for next feilvire location for a [O2/9O2]* SiCTRBSN laminate. 


that the assumption of regularly ^ced transverse cracks may be unrealistic (although the 
average crack location will always be the midpoint!). However, the global effects of the 
damage in these early stages are relatively minor. As damage progresses, and the crack 
drasity increases, the validity of the regular spacing assumption gains credence. 
Therefore, the assumption of a regular spacing throughout the damage development, which 
in turn allows the problem to be modeled via the analysis of a characteristic volume, would 
appear to be justified from an engineering standpoint. 
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Figure 53 : pdfs for next Mure location for a [0/904> Gi/Ep laminate. 


SUMMARY 

In this chapter, the method developed by Oh and Finney for modeling the characteris- 
tics of the location of faUnre has been revisited. Reworking the stress solution from 
Cluster 4 to accommodate the less general cross-ply system, an analytical egression was 
developed for the probability density fimction of the next failure location in a region of a 
cross-ply laminate bounded by two existing transverse matrix cracks. In exercising this 
model, it was found that, from a qualitative standpoint, the assumption of regularly spaced 
transverse matrix cracks is valid. 
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Summary 
and Conclusions 


Methcxls for probabilistically modeling the progressive failure of brittle matrix 
composite laminates were investigated in order to more accurately model the failiue 
process. Failure was approached from a non-interactive standpoint in which each of the 
in-plane stresses present in a composite layer were assumed to act independently of one 
another towards correspondingly separate and independent failure modes. 

A modified bundle theory which had been proposed in the literature to be a good model 
for the post-matrix-cracldng strength of a periodically cracked composite was examined and 
found to poorly predict the distribution in strength which has been commonly observed for 
these materials. Extensions to the model were introduced to allow for various material 
imperfections, in an effort to e^q)and the predicted distribution of str^gth. These remedies 
wm insufficient to solve the problem. Other possible extensions were indicated, but not 
investigated. 

A generalized shear-lag theory for analyzing the stress state of a transversely cracked 
layer having a general off-axis orientation and subjected to a state of full in-plane loading 
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was derived. This represents a major step in the development of shear-lag models, which 
prior to this had been limited in application to cross-ply laminates imder in-plane normal 
or shear loads (but not both). A method for utilizing this theory to model the effective 
bdiavior of a transversdy cracked layer within the laminate environment was established 
and a computer program was writtei to fully ocetcise die model. Results showing the 
characteristic stress-strain response, damage accumulation (in terms of the density of 
transverse matrix cracks) as a function of applied strain, and cumulative distribution curves 
for the probability of failure as a function of applied strain were presented for a variety of 
example problems. 

An analytical expression for the probability doisity fimction for the location of the next 
transverse crack occurrence within a crack bounded region was developed for a cross-ply 
laminate . The results of this investigation provided qualitative proof that the uniform crack 
spadng which had been assumed in the shear-lag analysis was a valid engineering 
assumption. 

This report has looked in detail at various ways by which macro-level behavior of 
composites may be extrapolated from micro-level analyses of rqiresentative damage 
r^mes. An approach of this manner should capture oiough of the actual physics of the 
itiiilure process to provide insight into the best ways to design with these matmial systems 
(this is sacrificed in a truly phenomenological model). At the same time the problem 
should be computationally simple enough to allow rq^plication of the model to practical 
RiTfd components (this not possible with a true micro-level, i.e. fiber/matrix, approach). 

In tile area of future work, experimental verification of the goieralized shear-lag model 
is needed. The results of Chapter 4 appear good qualitatively, but quantitative results are 
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necessary. Also, pertaining to the generalized shear-lag model, thermal stress effects need 
to be included. 
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Global Equilibrium 
Equations for 
Three Layer System 


The global equilibrium equations for a three-layer system were stated without proof in 
Chapt» 4 (equations 4.44, 4.45 and 4.46). In this Appoidix, the equilibrium conditions 
are derived in order to verify that the forms of these equations are in fact correct. 

The derivation begins by considering an undamaged three-layer laminate with 
dimensions and (see Figure A.l). This system is assumed to be subjected to 

an applied in-plane displacement field. The applied displacements results in a far-tield 
loading which is uniform (i.e. , the loads are independent of x, y and z ) and has a general 
in-plane nature (i.e., longitudinal, transverse and shear terms may be present). The far- 
fidd loads are specified in units of force per unit length. Due to the undamaged condition 
of the plate and St. Venant’s principle, the stress field which arises within the laminate as 
a result of the far-field loads is independent of x and y. However, because elastic 
prrq)erties can change from layer to layer, the stresses may vary from layer to layer as 
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Figure A.l: Three layer laminate with general in-plane loads applied. 

well, and are thus a function of z . Within each of these layers, in-plane loading terms are 
defined as the integrated stress through the layer thickness. 



(A.1) 

*1 



(A.2) 

*1 



(A.3) 


*1 


This integration has the effect of "smearing out" the z dependence of the stress, producing 
effective in-plane loads. 
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Consider the segment of the laminate resulting from slicing the laminate along the planp. 
x=x'. A free body diagram, showing the £ar-iidd loads on the free edges and the 
resultant in-plane loads on the exposed section for each of the layers, is presented in Figure 
A.2. Summing the forces in the x and y directions sq>arately, yields. 


EF, = 0 - -P,-L, - P^-x' * P^x- * 


P, - AT»> * AT- ^ W" 


(A.2) 

(A.3) 


£F, . 0 = -P, x- * P, x- - P^ L^ * (A.4) 


P = 

xy xy *^xy ^^xy 


(A.5) 


Similarly, sectioning the laminate along the plane y=y' (see Figure A.3) and summing 
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Figure A.3: Free body diagram for section along the plane y — y'- 
forces yields the equilibrium equation for the forces in the y direction. 

p, - w™ + w® + w® 

Equations (A.3), (A.5) and (A.7) comprise the global equilibrium equations, relating 
the fer-field loads to the resultant in-plane loads on each layer, for an undamaged laminate. 
The question which remains to be answered is whether or not these equations still apply 
for a laminate containing transverse matrix cracking. 

Due to the system being under a state of displacement control, the onset of damage will 
effect the magnitude of the far-field loads (as governed by the constitutive law), however 
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they are assumed to remain constant with respect to position (i.e., independent of x, y and 
z). Locally, the presence of damage within the laminate produces a stress field within 
each of the layers which is no longer independent of position. The assumption from 
Chapt^ 4 that damage effects are direction specific is continued here. For the case of 
transverse matrix cracks running in the x direction, damage is induced in the y direction. 
This implies that the stresses within the layers are functions of y and z, (i.e., 
of - ofiy-z)). 


"'fO’) - {<’!P(y.z)dz 

*1 


(A. 8) 
(A.9) 
(A. 10) 


and thus, the loading terms are functions of the y position at which they are evaluated. 

Returning to Figure A.3, a cut-away section of the laminate along the plane y=y' is 
examined. However, now the resultant in-plane loads are functions of y , Ny%y) and 
N^{y)y due to the transverse matrix cracks. Summing the forces on this segment in the 
X and y directions yields equilibrium equations for the transverse and shear forces. 


Py = w“(y') (a. ii) 

P:^ ‘ <’(>') (A- 12) 

With regard to the equilibrium in the x direction, a segment of finite length x' in thex 
direction, and differential length dy in the y direction (see Figure A.4) is examined. The 
far-field loads remain independent of position. The reaction forces on the faces y = y ' and 
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Figure A.4: Free body diagram for section of length x ' and width dy. 

y = y'+ dy are constant due to y being constant on these faces. The reaction forces on 
the exposed edge of this segment, jc =x\ vary with y . The average value of the in-plane 
load N^{y) over this region can be expressed in integral form as, 

Nf , T N^(y)dy (A. 13) 

dy 

Evaluating this integral using a first-order Taylor series expansion in y', it is found that 

^ (A. 14) 

If higher order terms are neglected, the load may be considered constant over the 
differential span. Similarly, the shear forces acting on this face can also be considered 

constant. 
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- <’(y') (A.15) 

Utilizing equations (A. 14) and (A. 15) and summing forces in the z direction yields the 
final equilibrium equation, 

P. = (A.16) 

Equations (A. 11), (A. 12) and (A.16) are the final form for the global equilibrium 
equations for a three layer system with damage present due to transverse matrix cracking. 
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Load Transformation 
of In-plane Loads 


Load terms are typicaUy thought of as vectors (or first order tensors) and are usuaUy 
specified in units of force [FORCE]. Accordingly, they transform via the first order tensor 

transformation equation. 


>1' 


1 

r 


>1^ 

p,. 

► [FORCE] = 

Oj.j 02*2 %'3 

< 

Pi 



/*3'1 ^'2 ^'3. 


Pi. 


[FORCE] 


(B.l) 


where the transformation matrix is comprised of the direction cosines a,.^- = cos{/ip^y). 
However, when referring to laminated plates, and plates in general, in-plane loads are 
usually specified in units of force per unit length of edge [FORCE/LENGTH]. These 
''loads’ are actually stress terms which have been integrated over a linear distance 
(specifically, the thickness of either the entire plate or an individual layer within the plate). 



(B.2) 

(B.3) 
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(B.4) 


Thus, what is typically ref«red to as the <^lied loading yector, , is actually 

a convenient grouping in column matrix fonn of elements which are actually members of 
a second order tensor, and must therefore be transformed accordingly. 

The transformation equation for a second order tensor is given in index notation as, 

(B.5) 

and in matrix notation as, 


[o'] = [A][o][Ar 

The transformation matrix, [A], is defined as 


(B.6) 


[A] = cos{A,^A.) H 


«i Pi 

”^2 P 2 (B.7) 

IJI3 «3 P3 

For an in-plane rotation as shown in Figure B. 1, evaluation of the direction cosines yields, 

s m = cos(6) 


«i = -«2 s « = sin(6) 

Pi = P2 = OT3 = /I3 = 0 

Pi = 1 

and results in the following form for the transformation matrix. 
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Figure B.l; In-plane transformation of axes. 


[^] 


m n 0 
-n m 0 
0 0 1 . 


(B.8) 


the matrix multipUcations of etjuation (B.6) for the case of an originaUy in-plane 


stress state, the transformed stress tensor is given by 


0^. T 


xY 




iVz Vx 


m n 0 
-n m 0 
0 0 1 


^xy ® 
^xy <»y 0 

0 0 0 


m n 0 
-n in 0 
0 0 1 


(m*o,+2mnr^+n^Oy) (-jnna,+(m^-n*)x^+inno,) 0 

{-mn0, + (m2-n2)x^+ntno,) (n^a^-2mnx^*m^Oy) 0 


(B.9) 


We are only concerned with those terms having non-zero values, i.e., o,,, Oy and 
These non-zero terms can be expressed in the following matrix equation. 
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nP’ iP 2mn 





tP nP -2mn 





-mn mn (jiP-tP) 




(B.IO) 


Similarly, the trHnsformEtion lelntion for the in-pl3ne lo^s is given by, 





fP 

2mn 


>x 

s 


nP 

-2mn 




-mn 

mn 






(B.ll) 


The reader should note that a matrix is simply a mathematical grouping device. It is not 
a tensor, and does not possess the physical characteristics of a toisor. As a result, the 
transformation of the in-plane loads must be carried out via equation (B.ll) and not 
through the first order tensor transformation equation. 
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Limitations of 
Generalized Shear-Lag 

Solution 


In Chapter 4, when equation (4.48) was reached in the derivation of the elasticity 
solution for the damaged region, the problem had been reduced to a second order ordinary 
differential equation containing two unknowns, and N^y (Note: this appendix reverts 
back to using the tilde notation to designate the homogenized system). This equation is 
repeated below in equation (C.l). 


2x7(1) 


d^N, 


y _ 


dy 


'2 


= ^ 




(C.1) 


In order to solve this equation, an additional independent equation was needed. For this 
purpose, the local equilibrium equation for layer 1 in the x' direction (equation (4.15)) 
was used, and combined with the shear-lag equations, and the thickness correction terms. 
This resulted in equation (4.53) (rewritten below as equation (C.2)) relating the second 
derivative of the shear load to the "difference in strain" term by the proportionality 
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constant i|r^, 



(C.2) 


where 


~ 2ih^+h ^)] 

By introducing global equilibrium considerations and diq>lacement continuity in the x' 
direction, this differential equation can be reformulated such that the strain terms are 
replaced by load terms (equation (4.55/C.4)) 


dy'^ 



(C.4) 


where o, ii and p have been defined in equations (4.49 - 4.51). 

For equation (C.4) to provide non-trivial information (i.e., for * 0), it is 

necessary to have a non-zero value for tlr^. Therefore, in order to determine the 
limitations of the solution given in Chapter 4 it must known under what conditions, if any, 

ti2 = 0. 

Looking at equation (C.3) the first obvious way in which a zero answer could be 
achieved is if ^ K ^ However, this is not a physical possibility for this problem 
set and can be disregarded. The second possibility would be if = = = 0 . 

To examine whether or not this occurrence is possible, the [H], [J], [IT] and [L] matrices 
are examined in greater detail. These matrices, originally defined in equations (3.28-31) 
of Chapter 3, are repeated below. 
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[L] > [[D]-[C][X]->[fl]]-‘ (C.5) 

[B] = [^]-U [A )->[«] [1][C][^)-' (C.6) 

[7] = -Mr‘[B][i) (C-7) 

m = -[i][C)(Xr‘ (C.8) 

where [A], [B], [C] and [D] are further defined in terms of the inverses of the layer 

stiffnesses. 

[-<] ■ - (C.9) 

[fl] = -y[<?i;‘ (C.10) 

[c] - -yraj’ (C.11) 

[D] ' (C-12) 


The stiffness matrices appearing in the above equations relate the out-of-plane shear stress 


to strain for each of the layers, i.e., 



i = i,2,3 


(C.13) 


with the primes signifying that the terms have been transformed to the x'-y' coordinate 
system. 

If each of the stiffness matrices involved in the above computations were to become 
diagonal matrices, i.e. , 


[Qlr 


Q'. 


55 


0 Q\ 


44 


i = 1,2.3 


(C.14) 


then in all subsequent calculations , the off- 
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diagonal terms would also be zero. Therefore, in order to determine under what conditions 

= 0, it is necessary to determine the circumstances leading to the 
terms taking on zero values for all three layers at the same time. 

The transformed coordinates x'-y' are aligned with the material coordinates of layer 
1 , therefore [<?']i = [<?]j • Since the material is orthotropic, = 0 wiU always hold 
true. 

Layers 2 and 3 are homogenized layers. Their elastic properties come from averaging 
the elastic properties of all the layers (except for layer i , the layer being analyzed for 
transverse matrix cracking) from the original system. 


[Q'h = [Q'h 





(C.15) 


The stiffness term <?;, is an odd function of the angle of orientation, that is to say 
Q^+ 0 ) = -Q^-6 ) . Therefore, if in the transformed coordinate system the orientations 
of the undamaged layers are balanced with respect to layer i , the summation of the Q'^s 
for the individual layers will cancel each other out and sum to zero. Thus, if the laminate 
configuration is balanced with respect to the damaged layer, the off diagonal stiffness term 
will equal zero. As a result, will tate on a zero value for this case. 

Additional possibilities for which achieve a zero value are if either 

= Fu'M * ® waetoo«”- Using the 
symbolic manipulator MACSYMA*, these expressions can be evaluated in terms of the 
inverted stiffness matrices for each layer. The resulting expressions are extremely 
complicated in an algebraic sense. The complexity of these expressions makes it 
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impossible to say without reservation whether either of the above conditions could occur, 

although it can be stated that it appears unlikely. 

Thus, the Umitation to the solution of Chapter 4 is that it is not appUcable to laminates 
whose configuration is such that the undamaged layers are balanced with respect to the 

damaged layer. 
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The PFRAC Algorithm 
for Laminate Analysis 


The computer code PFRAC, Progressive Failure and Rdiability Analysis of 
Composites, which is used in this rqx}rt models the progressive failure of composite 
laminates under strain-controlled conditions. The material strengths are considered to be 
random variables of known distribution, and multiple modes of failure (i.e., longitudinal, 
transverse, and shear (see Chapter 1)) are considered within each layer. As individual 
modes foil, corresponding stiffness reductions are made via the shear-lag methods of 
Ch^ters 3 and 4, as well as traditional ply drop-off (or in this case mode drop-off ) 
techmques. Monte Carlo methods are employed to determine the cumulative distribution 
curves for the laminate’s probability of failure as a function of the applied strain. This 
Appendix provides an overview of theory and logic behind the PFRAC code. 


SIMULATING THE STRAIN CONTROLLED TEST 


The PFRAC code simulates a strain-controlled environment. The laminated plate in- 
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plane constitutive relationship is given by (Jones, 1975) 



A' 

H' 



(D.l) 


with €* being the mid-plane strains, k being the plate curvatures, and JV and Af being the 
in-plane loads and moments, respectively. The prescribed conditions on the system to be 
inodei«t are: a known value for the applied axial strain, stress-free states for the 

in-plane transv»se and shear loads = 0 and = 0), and moment-free in-plane 

conditions (Af, = 0). Substituting these known conditions into equation (D.l), 

the unknown strains and loads are given by the following expressions. 



(D.2a) 

K ^ 

II 

(D.2b) 

II 

(D.2c) 

can be determined as 

a function of the applied 


longitudinal strain. 

With the full in-plane strain vector known, the stresses within each layer can be 
detramined 


' * 


€x 


• -[Q]r 



i 

.Yx,. 


(D.3) 


In equation (D.3), [<?]j represents the stiffness matrix of layer i transformed to the global 

coordinate system. These individual layer stresses can, in turn, be transformed back to the 
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material coordinate systems for a specified layCT, and these stresses can be compared to 
the corresponding material strengths, X, in order to evaluate failure (i.e., / i 1). 



(D.4a,b,c) 


Based on the stress solution for a very small applied strain, the program can determine 
which failure mode will be the first to fail and st^ to the applied strain level correspond- 
ing to its failure. This is possible due to die linear bdiavior of the undamaged composite. 
The strain is increased incrementally from this point forward. 

Once a failure state has been attained in an individual mode, steps must be taken to 
reduce the stiffness of that mode in an appropriate manner. Failure of modes 1 and 6 are 
one time occurrences; that is, once they fail, their corresponding stiffnesses (£j for mode 
1 and Gj 2 for mode 6) are reduced and they are dropped from future computations all 
together. (Note, the stiffnesses are reduced by a fector of 10“^ rather than being set equal 
to zero. This is done to avoid encountering singularity problems in succeeding computa- 
tions.) A mode 2 failure is modeled using the shear-lag theory which allows for the 
damage evolution to be included rather than having complete damage being imposed 
immediately as is done with the drop-off techniques. 

The shear-lag model determines the effective elastic properties of the damaged layer 
based on the homogenized elastic properties of the rest of the laminat e and the present 
value of the in-plane strain vector { e,, }’’. However, referring to equations (D.2), 

it is apparent that the in-plane strains are calculated as a function of the laminate stiffness. 
The new predicted values for the effective stiffness of the damaged layer will influence the 
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updated values for the in-plane strains. Thus, an iterative process is required. The 
following cost function is set up, 

where e is the previous strain value, and c* is the updated strain value. Note that e, does 
not zppeax in the cost function because its value is known and remains fixed until 
incremented to the next strain level. Theoretically, the correct solution has been reached 
when the cost function equals zero. Within PFRAC, the convergence criteria has been set 
to C < 1 X 10"* . Experience indicates convergence to this value happens quickly , and this 
particular value returns relatively good accuracy for the corresponding strain values. 

When the interfacial shear stress reaches its Umit, and no increase in load transfer is 
possible, then this mode is considered to be finaUy failed. However, while no additional 
load can be sustained by the layer, it is assumed that the layer continues to carry the load 
that it presently bears. This assumption is accounted for by storing the current layer load 
in a buffer and reducing the stiffnesses of the layer by a factor of 10"’ , so as to effectively 
remove the layer from all future computations. The loads contained in the buffer are then 
^jpended to the computed laminate loads when evaluating the total load supported by the 

laminate. 

After all the necessary stiffness reductions have been completed, the [A], [B] and[D] 
matrices of the laminate constitutive equation are reformulated, and subsequently the stress 
solution is updated. If additional failures are encountered under the updated stress state 
(as in a "domino effect"), the appropriate stiffness reductions are made in the manner just 
presented. If no additional failures have arisen, the applied strain is incremented, and the 
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stress state is evaluated for the new appUed strain. The whole process continues until total 
laminate failure is reached. There are two criteria by which total laminate failure is 
judged. The first criterion accounts for a sudden drop in the load bearing capacity of the 
laminate. If, for a given applied strain, internal fidlures occur resulting in a decrease in 
the load by 50% or greater, criterion 1 has beai met. A gradual decrease (i.e., over a 
range of strain values) of 50% is allowable, only a suddra decrease foils the laminate. The 
second criterion deals with the positive definiteness of the material stiffness matrix. The 
stiffness matrix for any material is required to be positive definite in order to ensure that 
the strain-energy density function be closed and convex. Thus if the effective stiffness 
matrix calculated to satisfy the shear-lag model no longer abides by this requirement, the 
laminate is considered failed. 

PROBABILISTIC CONSIDERATION 


It is assumed that the material strengths, X. (i = 1, 2 and 6) corresponding to each of 
the individual failure modes are characterized by a Weibull distribution, i.e.. 



(D.6) 


where the distribution parameters ( characterizing each failure mode can be obtained 
from strength data. Due to the random nature of the strengths, the reliability for an 
individual hypothetical sample will have a value between zero and one inclusive such that 
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Fx^ = \-u (D.7) 

where u is a uniform random variable contained on the interval [0,1]. 

Substituting equation (D.7) into equation (D.6) and inverting the distribution function 
yields the following realization of the random strength variable, 

X. = p^[-ln(l-ii)p (D.8) 

Here represents a sample value of X,. . By utilizing a uniform random number generator 
to obtain values of u , realizations of the strength can be obtained. Within the Monte Carlo 
procedure, this is done many times in order to achieve a sample population for the 
material . This sample population is analyzed for failure in the manner described in the 
previous section, and the applied strain present at failure is saved. The recorded failure 
strains are then ranked and assigned a corresponding probability of failure, p , using a 
median ranking 

Pj = ; j = 1.2.3....,n (D.9) 

^ n+OA 

wh^ n is the total number of samples. 
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